Orbit Determination with the two-body 

Integrals. II 

Giovanni F. Gronchi, Davide Farnocchia, Linda Dimare 

Dipartimento di Matematica, Universita di Pisa 
gronchiOdm.unipi . it 
f arnocchiaSmail . dm . unipi . it 
dimareOmail . dm . unipi . it 

Abstract 

The first integrals of the Kepler problem are used to compute pre- 
liminary orbits starting from two short observed arcs of a celestial 
body, which may be obtained either by optical or radar observations. 
We write polynomial equations for this problem, that we can solve 
using the powerful tools of computational Algebra. An algorithm to 
decide if the linkage of two short arcs is successful, i.e. if they belong 
to the same observed body, is proposed and tested numerically. In 
this paper we continue the research started in |6j, where the angular 
momentum and the energy integrals were used. A suitable component 
of the Laplace-Lenz vector in place of the energy turns out to be con- 
venient, in fact the degree of the resulting system is reduced to less 
than half. 

1 Introduction 



We present a new method, based on the first integrals of the Kepler problem, 
to compute a finite set of preliminary orbits of a celestial body from two short 
arcs of observations. We assume that the body moves on a Keplerian orbit 
with a known center of attraction ojll and is observed from a point P, whose 

-'^For asteroid orbits O corresponds to the center of the Sun, for space debris O is the 
center of the Earth 
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motion is a known function of time. We deal with two different kinds of 
observations, optical and radar, and we make use of the related attributables 
(see ig, [T2])0 

In [6] the angular momentum and the energy integrals are used to solve 
the linkage problem for solar system bodies. This means to identify two 
attributables as related to the same observed object by computing (at least) 
one reliable orbit from the observations of both attributables. The equations 
of the problem are written in a polynomial form and the total degree of the 
system is 48. The use of these integrals for the linkage problem has been first 
proposed in [IQj, but without fully exploiting the algebraic character of the 
problem. 

The algorithm presented in [B] has been used in [3] for the problem of 
correlatior^ of space debris: here the authors have extended the method 
including the oblateness effect of the Earth. 

In this paper we propose different equations for the same problem: in 
particular we use a suitable projection of the Laplace-Lenz vector in place 
of the energy. The advantage of this approach is that there are several 
cancellations and the total degree is 20. 

The same equations can be written using different data, simply consider- 
ing other quantities as unknowns: in Section |4] we deal with the case of an 
optical and a radar attributable. This case is peculiar because we end up 
with a univariate polynomial of degree 4. Thus this problem admits explicit 
solutions. 

In both cases the solutions must fulfill compatibility conditions (as also 
shown in [6]), taking into account the other integrals of Kepler's problem. To 
select the solutions we propose a different strategy, based on the attribution 
algorithm of a very short arc to a known orbit, see [H], [Z]- 

The structure of the paper is the following. After introducing some def- 
initions in Section [21 we study the linkage of two optical attributables in 
Section [3l while in Section H] we consider the same problem with one optical 
and one radar attributable. The degenerate cases are shown in Section |5l 
Sections E] and [7] are devoted to explain the computation of the covariance 
matrix for each orbit and the selection of the solutions. We conclude with a 
numerical test in Section ID 

^The two different attributabies can be obtained from observations made from different 
stations 

■^that is tlic linkage problem, in tlie context of space debris 
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2 Preliminaries 



Let us fix an inertial reference frame, with the origin at the center of attrac- 
tion O. The position q and velocity q of the observer are known functions 
of time. We describe the position of the observed body as the vectorial sum 

r = q + pe^ , (1) 

with p the topocentric distance and e'' the line of sight unit vector. We 
choose spherical coordinates (a,(5, p) G [— 7r,7r) x (— 7r/2,7r/2) x R+, so that 

e'' = (cos 6 cos a, cos 6 sin a, sin 6) . 

A typical choice for a, 6 is right ascension and declination. Then we can 
write the velocity vector 

r = q + pe^ + p(d cos (5e" + (5e^) , p,a,6 eR, p eR^^ , (2) 

where p, pa cos 6, p6 are the components of the velocity, relative to the 
observer in P, in the (positively oriented) orthonormal basis {e'',e",e^}, 
with 

We recall the definitions of optical and radar attributables. From a short 
arc of optical observations of a moving body (tj,aj,5j) with i = l...m, 
m > 2, it is possible to compute an optical attributable 



Aopt = ia,5,a,S) e [-n,n) x {-n/2,n/2) x 



representing the angular position and velocity of the body at a mean time t 
(see |8],[6]). In this case the radial distance and velocity p, p are completely 
undetermined and are the missing quantities to define an orbit for the body. 
From a set of radar observations of a moving body {ti,ai,6i, pi), with i = 
1 . . . m, m > 2, it is possible to compute a radar attributable, i.e. a vector 

Arad = (a,5,P,p) e [-7r,7r) X (-7r/2,7r/2) x M+ x M, 

at time i (see Here a, 6 are the unknowns needed to define an orbit. 

We call attributable coordinates the vector {a,6,a,6, p, p) representing 
the position and velocity of the body as seen from the observer at time t. 
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3 Linking two optical attributables 



Given two optical attributables Ai, A2 at epochs ti, ^2, we assume they belong 
to the same observed body and write 4 scalar algebraic equations for the 
topocentric distances pi, p2 and the radial velocities pi, p2 at the two epochs. 
We use some of the algebraic integrals of the Kepler problem, i.e. the angular 
momentum c, and the Laplace-Lenz vector L. The expressions of these 
integrals as functions of the topocentric distance and radial velocity p, p are 
given below. 

Angular momentum: 

c(p, p) = r X r = Dp + Ep^ + Fp + G , 

where 

D = q X e'' , 

E = d cos 6eP x e" + Se^ x = a cos ^e*^ — 6e°' , 
F = a cos 5q X e" + (5q X e"^ + e'' X q , 
G = q X q . 

LaPLACE-LeNZ'S VECTOR: 

lxL{p, p) = r X c — p— = ^|r|^ — ■pj)'^ ~ ' ■> 

where p is a positive constaniQ and 
|r| = (p2 + |q|2 + 2pq. 6^1/2^ 

|r|2 = p2^(d2cos2 5 + 52)p2^2q■e^p + 2q•(acos(5e° + 5e'^)p+|q|^ 
r ■ r = pp + q • e'^p + (q ■ e'' + q ■ e^d cos 5 + q ■ e^5)p + q ■ q . 

Remark 1. If O corresponds to the center of the Sun, then we use interpo- 
lated values for q, q, as suggested by Poincare If O corresponds to the 
center of the Earth we do not apply this method. 

These dynamical quantities give 6 scalar integrals of the motions: only 5 are 
mutually independent, in fact we have L ■ c = 0. Since we have 4 unknowns, 
generically we only need 4 scalar conservation laws to define a finite number 
of solutions. We select the conservation of the angular momentum vector 

^/i = GrriQ if we deal with objects orbiting around the Sun; /i = Gm^ for satelhtes of 
the Earth. 
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and of a particular component of the Laplace-Lenz vector. The choice of the 
latter integral presents a substantial advantage with respect to the use of 
the energy, as in [6] : the difference between the two choices will be discussed 
later. 



3.1 The polynomial equations 

We use the notation above, with index 1 or 2 referring to the epoch. If Ai, 
A2 correspond to the same observed object, then the angular momentum 
vectors at the two epochs must coincide: 

Cl(pi,Pl) = C2(P2,P2) • (3) 

Equation ([3]) can be written as 

Dipi - D2P2 = J(pi, P2) , (4) 

where 

J(pi, P2) = E2p^ - Eip2 + F2P2 - Fipi + G2 - Gi . 
Following [6j we eliminate the variables pi,p2 and obtain the equation 

Di xD2- J(pi,p2) = . (5) 

We can write the left-hand side of as 

q{pi, P2) = q2opl + Qiopi + %2pl + goiP2 + goo , (6) 

with 

g2o = -El ■ Di X D2 , go2 = E2 ■ Di x D2 , 

gio = -Fi • Di X D2 , goi = F2 ■ Di x D2 , 

goo = (G2 - Gi) • Di X D2 . 

The radial velocities are given by 

. (J X D2) • (Di X D2) . (J X Di) • (Di X D2) 

Pi(pi,P2) = |D^^D^|2 ' P2(pi,P2) = |D^^D^|2 • (7) 

Also the Laplace-Lenz vectors at the two epochs must coincide. We equate 
the projection of both vectors along v = 63 x q2: 

Li(pi,pi) ■ V = L2(P2,P2) ■ V . (8) 
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Actually the projection of L2 along v is particularly simple: 



/iL2 ■ V 



1-2 ■ r2)(r2 ■ v) 



thus dS]) becomes 



( 




(ri ■ v) - (ri ■ ri)(ri ■ v) = -(1-2 ■ r2)(r2 ■ v) . (9) 



After substituting ([7]), this is an algebraic equation in pi,p2- Rearranging 
the terms in ([9]) and squaring we obtain 



p{p^,P2) =^ ^'(ri • v)2 - { [Iripn - (ri • ri)ri + (r2 • r2)r2] • v}' = . (10) 



This is a polynomial equation of degree 10 in pi, p2'- in fact, the projection 



does not depend on p2 and, in the difference |ripri — (ri ■ ri)ri, the second 
degree term in pi (i.e. pfpieD cancels out. 

Therefore, to solve the linkage problem, we can consider the polynomial 
system 



with total degree 20. This shows the advantage of this method compared 
with the one in [H], which gives total degree 48. 

3.2 Computation of the solutions 

To compute the solutions of f fT2|) we define an algorithm similar to the one in 
0' 0- -^y grouping the monomials with the same power of p2 we write 



1-2 • V = q2 • (p2(a2 cos ^262 + ^262) + q2) x = 

= P2(-a2 cos52q2 ■ - '^2q2 ■ ) + ■ q2 X q2 



(11) 




(12) 
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where 



(13) 



j=0 




10 for j = 

10 - (j + 1) for j = 2A; - 1 with k>l 
10 — j for j = 2k with k > 1 
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and 



q{pi,p2) = &2 P2 + ^1 P2 + boipi) 



(14) 



for some univariate polynomial coefficients 0^,60 and constants &i,62- 
We consider the resultant Res{pi) of p, q with respect to P2: it is generically 
a degree 20 polynomial defined as the determinant of the 10 x 10 Sylvester 
matrix 



f as 



S(pi 
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bo 
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bi 
bo 
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bi 



V ao 
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bo bi 
60 / 



(15) 



The positive real roots of Res{pi) are the only possible values of pi for a 
solution (pi,p2) of f|T2|) . 

Remark 2. T/ie resultant ofp,q with respect to pi leads to compute deter- 
minants of 12 X 12 matrices, thus the elimination of p2 is more convenient. 
On the other hand, if we project the Laplace-Lenz vectors on x qi, it is 
better to eliminate pi . 

Following P, we compute the coefficients of Res{pi) by an evaluation- interpo- 
lation method based on the FFT, and then the roots pi{k) of Res{pi) by the 
algorithm described in |lj. The computation of the preliminary orbits is 
concluded as follows: 

1) solve the equation q{pi{k), P2) = 0; 

2) discard spurious solutions, that is pairs {pi, P2) solving ffTOj) but not 



3) compute the values of pi{k),p2{k) by ([7]); 

4) write the corresponding orbital elements. The related epochs are ti = 
ii — Pi{k)/c, i = 1,2 where c is the velocity of light (aberration correc- 
tion). 
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4 Linking radar and optical attributables 



Assume we have a radar attributable Arad = {c(,S,p,p) at epoch t. We 
introduce the variables 

^ — pa cos 5 , C — 

so that 

r ^ + Ce^ + (pe^ + q) , 

r • r = q • + q • e'^C + (pe'' + q) • r . 
The angular momentum as a function of ^, C is 

CradiCC) = AC + BC + C, (16) 

where 

A = rxe", B = rxe'', C = rxq + pqxe''. 

Suppose wc have a radar attributable Arad at time ti and an optical 
attributable Aopt at time t2- Equating the angular momentum vectors c^ad 
and Copt at the two epochs we obtain a polynomial system of 3 equations in 
the 4 unknowns CiXi^ P2, P2' 

AiCi + BiCi + Ci = D2P2 + E2PI + F2P2 + G2 . (17) 

The system is linear in ^1, Ci, P2- By solving for these variables we obtain 

6(P2) = X2PI + X1P2 + X0 

Cl(p2) = Z2P^ + ZiP2 + Zo , (18) 

Mp^) = R2pi + R1P2 + Ro 

where 

X2 = 7E2-Bi XD2, Xi =7F2 -Bi X D2, Xq = 7 (G2 - Ci) • Bi x D2 , 

Z2 = -7 E2 • Ai X D2 , Zi = -7 F2 • Ai X D2 , Zo = -7 (G2 - Ci) • Ai X D2 , 

R2 = -7E2-Ai xBi, Ri = -7F2-Ai xBi, Rq = -7 (G2 - Ci) • Ai x Bi , 

and 7 = l/(Ai • Bi x D2). 
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Equating the expressions of the Laplace-Lenz vectors at the two epochs, and 
projecting along v = eg x q2, yields 

/i[Lrad(6, Cl) - Lopt(P2, P2)] ■ V = 

= (^|riP - i^jri - (ri • ri)ri ■ v + (1-2 ■ r2)(r2 ■ v) = . 
The term 

1-2 ■ V = p2(-a2 COS 62^2 ■ el + 62^2 ■ + 63 ■ q2 x q2 

does not depend on p2 and is linear in p2 (cfr. with (TTTl)). Thus, after substi- 
tuting p2 = P2(P2), 6 = 6(P2),Ci = Ci(P2) from ([IHD, the terms (r2-r2)(r2-v) 
and [|ripri — (ri ■ ri)ri] ■ v are polynomials of degree 4 in p2- We obtain 
a univariate polynomial equation with degree 4 in p2; which admits explicit 
solutions. For each positive root p2{k) we can compute orbital elements at 
epochs h = ti- Pi/c, ^2 = h - P2{k)/c using (IH]). 

5 Degenerate cases 

We list the cases that make the equations of the linkage degenerate. 

Optical case 
The quadratic form (EI) is completely degenerate if 

El ■ Di X D2 = E2 ■ Di X D2 = . 

For a discussion on the geometric meaning of these conditions see [6]. Another 
degenerate case occurs if 63 x q2 = 0. In the case of space debris this 
corresponds to a zenith observation. 

Radar-Optical case 
System (IT7|) degenerates if 

Ai X Bi • D2 = r^(ri ■ D2) = . 

This occurs when ri ■ e^* = 0, or ri x r2 = 0, or when 63 is in the orbital plane 
(orthogonal to ri x r2). Another degeneration occurs when 63 x q2 = 0, as 
in the optical case. 
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6 Covariance of the solutions 



Let A = (^1,^2) be the vector of two optical attributables and Fa its co- 
variance matrix. For each solution Y = (pi, pi, p2, P2) of the linkage problem 

r c,(p„pO = C2(P2,P2) 

[ Li(pi,pi) • V = L2(P2,P2) ■ V 

we can compute the Cartesian coordinates Ecm--, Sml at epochs ti,t2, and 
their covariance matrices T^r, F^ar- We introduce the following notation: 

1) Ecar- = {Scar, Seal) is the 2 cpochs Cartcsiau coordinates vector; 

2) E,,i = wher^S 

S^itt ~ i'^i} ^11 Pii Pi) 1 i = 1, 2 . 



Define the map ^ : 



)12 



— )■ 



by 



Cl - C2 

p(Li - L2) ■ w 



w = r2 X q2 



Eatt — Ecar by ([T]) , ([2]) for both epochs, and consider 



Moreover, define T^t, „ „ , „ 

the map $ = * o T^^f . Then $ = is equivalent to ([T^dFI 
The covariance matrix of the Cartesian coordinates at epoch ti is 



r(i) 



(1) 

car , 



dA 



(1) 

car 



with 



(1) 

car 



(1) f)cW 
car ^<--att 



OA 



From the implicit function theorem 



dY 

dX' 



A) 



dS. 



(1) 

att 



dA 



h O4 
d{pi,Pi) 
dA 



E 



att) 



dA^ 



E, 



att) 



^If we use interpolated values for q, q, as suggested in |9] , then £^l{ are not the at- 

(i) 

tributable coordinates corresponding to Scar, i— 1,2. 

^We use w instead of v to obtain simpler expressions for the derivatives of 
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where 
dY 



°latt 



dE 



car 
att 



dY ' 



dA 



„ '-f-car 
°latt 



car 
att 



dA 



The matrices and are respectively made by columns 5,6,11,12 and 



by columns 1,2,3,4,7,8,9,10 of 
For a given vector u e define the hat map 



dEatt ■ 



3 {Ui,U2,Us) = U 1-^ U 



def 






-Us 


U2 









-Ml 


e so(3) 


-U2 


Ui 








Then we have, using 
(9* 



dEr 



-u, 

-ri Yi 1-2 -r^ 

d^c dAc dAc dAc 

dri dvi dr2 dr2 



where 

dAc _ 

dri 
dAc 

dri 
dAc 

dT2 

dAc 

dY2 



T , (ri • w) J, J, 
— ri - (ri • w)ri , 



ki / l^i 



= 2(ri • w)rf - (ri ■ w)rf - (ri • ri)w^ , 



|ri| 



[q2 X ri]'^' - (ri • ri)[q2 x ri]'^' + (h • w)r2 + • r2)[q2 x Y2Y' , 



(1-2 • w)r2 + (r2 • r2)w . 



In the case of one radar and one optical attributable the covariance of the 
solutions can be computed in a similar way, with the following differences: 

1) the vector of the attributables is A = {Arad, ^opt); 

2) the vector of unknowns is Y = (di, 5i, p2, P2), 

3) the matrix of the derivatives of S^u with respect to A is 

I2 O2 O2 O2 
d{ai,6i) 



(1) 

att 



dA 



dA 

O2 h O2 O2 
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4) The matrices and are respectively made by columns 3,4,11,12 
and by columns 1,2,5,6,7,8,9,10 of 

7 Selecting solutions 

The solutions of f[T^ are defined by using only four conservation laws. Thus 
Sell, Seal may not correspond to the same orbit. We select the solutions of 
the linkage problem by means of the attribution algorithm [8] , [7] . Here we 
recall briefly the procedure. 

Let £i be a set of orbital elements for the observed body at time ti, with 
6x6 covariance matrix Fi. We can propagate the orbit with covariance to 
the epoch t2 of an attributable A2, with a given 4x4 covariance matrix F^j, 
by the formula 



1 



dSi 



z \-\T 



where ^{Si,t) is the integral flow of the Kepler problem. Then we can 
extract a predicted attributable Ap, at time ^2, with covariance matrix F_4p. 
Let C^^ = (F^J-i and C^, = (F^J-^ We define 

Co = + C*^2 5 To = Cq ^ 
The identification penalty is given by 

= (A2 - Ap) ■ [Ca, - c^,roC^j(A - Ap) . 

If the value of X4 is within a fixed threshold, we can accept the orbit Si. 

Remark 3. To select solutions we could also use compatibility conditions, 
as in fSjl- In this case the conditions could be 

(Li-L2)-e^ = 0, ^i-^2 = ni(ti-t2) . (20) 



8 A test case 

We present the results of a test of the method explained in Section [3] with 
the asteroid (99942) Apophis. We take two sets of 13 and 12 observations 
respectively with mean epochs ti = 53175.59, t2 = 53357.45. After removing 
duplicate and spurious solutions we obtain 
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Pi 


P2 


1 


0.78987 


0.04345 


2 


1.13777 


0.09569 



Table 1: Solutions of the system ^ for ((99942)). 

The two solutions gives respectively X4(l) = 3230925.94, X4(2) = 2.29, there- 
fore we select the second one, with Keplerian elements (distances in AU, 
angles in degrees) 

a = 0.9230 , e = 0.189 ,/ = 3.287 , n = 204.912 , to = 124.778 , i = 249.003 

at epoch ti = 53175.59. We can compare the results with the known orbit 
propagated at epoch ti. 

a = 0.9219 , e = 0.191 , / = 3.333 , n = 204.575 , to = 126.176 , i = 247.500 

In Figure [H for the test case of (99942) Apophis, we show the intersections 
of the curves defined in this paper compared with the ones obtained by 
the conservation of the energy. In the four pictures the hyphened curve 
corresponds to equation ([5]). We also draw the curve defined by ([9]) on top 
left, and the one by ( ITOl) . in polynomial form, on top right. The conservation 
of the energy defines the curve drawn on bottom left, its polynomial form 
(obtained by rearranging terms and squaring twice) defines the one on bottom 
right. The orbit determination method introduced in this paper, searching 
for the intersections shown on top right, is clearly convenient with respect to 
the method investigated in [0], related to figure on bottom right. 

9 Acknowledgments 

We wish to thank A. Milani for his useful suggestions during the development 
of this work. 



13 



Figure 1: For the test case of (99942) Apophis, this figure shows the advantage 
of using equation ([5]) instead of the conservation of the energy £. Top left: c, 
L ■ V integrals. Top right: c, L • v integrals, polynomial form. Bottom left: c, £ 
integrals. Bottom right: c, 8 integrals, polynomial form. 
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